Evaluation of ground state entanglement in spin systems with the random phase 

approximation 
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We discuss a general treatment based on the mean field plus random phase approximation (RPA) 
for the evaluation of subsystem entropies and negativities in ground states of spin systems. The 
approach leads to a tractable general method, becoming straightforward in translationally invariant 
arrays. The method is examined in arrays of arbitrary spin with XY Z couplings of general range in 
a uniform transverse field, where the RPA around both the normal and parity breaking mean field 
state, together with parity restoration effects, are discussed in detail. In the case of a uniformly 
connected XY Z array of arbitrary size, the method is shown to provide simple analytic expressions 
for the entanglement entropy of any global bipartition, as well as for the negativity between any 
two subsystems, which become exact for large spin. The limit case of a spin s pair is also discussed. 

PACS numbers: 03.67.Mn, 03.65.Ud, 75.10.Jm 



I. INTRODUCTION 

The study of entanglement constitutes one of the most 
active and challenging research areas, being of central in- 
terest in the fields of quantum information [l[ and many- 
body physics 0. The concept of entanglement has pro- 
vided a new perspective for analyzing quantum correla- 
tions and quantum critical phenomena in many particle 
systems, leading to fundamental results and new insights 
in the field 0t3- Nonetheless, the evaluation of entan- 
glement in general strongly interacting many-body sys- 
tems remains a difficult task, particularly in systems with 
long range interactions, high connectivity and large di- 
mensionality, where usual treatments such as Quantum 
MontcCarlo @, DMRG or matrix product states [|[ 
become more involved or difficult to implement. In pre- 
vious works 0, [n| we have applied a general mean field 
plus RPA treatment to the evaluation of pairwise entan- 
glement (i.e., that between two elementary components) 
in spin systems at zero and finite temperature. The ap- 
proach was able to capture the main features of the en- 
tanglement between two spins in arrays with XY and 
XYZ couplings of different ranges, including the predic- 
tion of full range pairwise entanglement in the vicinity 
of the factorizing field [Tol - fl2l |. The accuracy of the ap- 
proach was shown to increase with the interaction range 
or connectivity. 

The aim of the present work is to examine the ca- 
pability of the previous method for predicting, in the 
ground state of spin systems, the entanglement proper- 
ties of general subsystems. We will focus on the entan- 
glement entropy of arbitrary bipartitions of the whole 
system, as well as on the negativity between any two sub- 
systems, not necessarily complementary, where the rest 
of the spins play the role of an environment and entangle- 
ment can no longer be measured through the subsystem 
entropy. Other measures, like the negativity (an entan- 
glement monotone computable for general mixed states 
US EH) have to be employed. This type of entangle- 
ment has recently received special attention [r^ - UtT ] since 



its behavior can differ from that of global bipartitions. 
We will show that the present approximation provides a 
general tractable scheme for evaluating these quantities, 
becoming analytic in translationally invariant systems. 

In Section II we present the general RPA formalism, 
describing the RPA spin state, the associated bosonic es- 
timation of subsystem entropies and negativities, the im- 
plementation in translationally invariant systems and the 
application to a general spin s array with XYZ couplings 
of arbitrary range in a transverse magnetic field. Symme- 
try restoration effects in the case of parity-breaking mean 
fields are also discussed. As illustration, we derive in sec. 
IIII results for a spin s pair and for a fully connected 
finite spin s array, where RPA is able to provide sim- 
ple full analytic expressions for subsystem entropies and 
negativities, which represent the exact large spin limit at 
any fixed size. Appendix lAl discusses the equivalence be- 
tween the spin and the bosonic RPA treatments, whereas 
appendix IBI contains details of the analytic results of sec. 
Conclusions are drawn in IV. 



II. FORMALISM 

A. RPA for spin systems at T = 

We will consider a general finite system of spins = 
{six, Siy, Si Z ), connected through general quadratic cou- 
plings and immersed in a magnetic field, not necessarily 
uniform. The corresponding Hamiltonian is 



H = B^i 



I . it 



(1) 



where fi = x,y, z and B 1 ^ are the field components at 
site i. Ising, XY, XYZ (J^ = S^Jjf) as well as 
Dzyaloshinskii-Moriya ( J l ^-'J V = — J W J>) couplings of ar- 
bitrary range are particular cases of Eq. ([1]) . 

The first step in the RPA [lj| is to determine the mean 
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field ground state, i.e., the separable state 
|0> = ®?=i|0i) = |0i...0„) 
with the lowest energy (H)o = (Q\H\Q), given by 

(H) =J2 Bitt ( s ^)o-l E ^"{'Msjvh (2) 

where (s,)o = (0,|si|0j). Each local state |0j) can be 
determined self-consistently as the lowest eigenstate of 
the local mean field Hamiltonian 

7 _ 9 ( H )o _ \i „ fo\ 

„ 9(si M )o 

being then the state with maximum spin Si directed along 
—A' (a local coherent state). This leads to the self- 
consistent equations 

A v = B v _ J^ jv (s jv ) , ( Sl ) = -SiX/X , (4) 

where A' — |A'|. Eq. (|3|) can be solved iteratively start- 
ing from an initial guess for |0,) or Aj, although other 
procedures (like the gradient method) can be employed. 
Eq. © becomes then (H) = \ J2i( xi + ' ( s i)o- 

Since the form ([!]) is valid for any choice of the local 
axes, it is now convenient to choose Zi along \ l , such 
that (si M ) = SiS^x and A v = X l S^ z , with X 1 > 0. The 
second step in the RPA is the approximate bosonization 

s i+ \/2slb\ , s,_ -> \f2sih , s iz -s l + bj^ , (5) 

where Sj± = si x ±isi V and 6i, &j are considered standard 
boson operators ([bi, &t] = 5^ 6j] = \b\, &t] = 0), with 
|0) — > I Oft) their vacuum. This bosonization is in agree- 
ment with that implied by the path integral formalism 
of [H, [ltj for T — s- 0, and preserves two of the exact spin 
commutators exactly ([s*,sf] = zLSijsf), the remaining 

one preserved as vacuum average (([s^,Sj])o = 2sj5jj). 
It coincides with the Holstein-Primakoff and other exact 
bosonizations [l8l - [2"lj up to zeroth order in s^ 1 . 

The third step is to replace Eq. J5]) in the original 
Hamiltonian ([I]), neglecting all cubic and quartic terms 
in bi, b\. This leads to the quadratic boson Hamiltonian 



linear terms in bi, b\ appear in H b , reflecting the stabil- 
ity of the mean field state |0) with respect to one site 
excitations. 

The last step is the diagonalization of the bosonic 
quadratic form ([6]), which is always possible when the 
hermitian matrix H. in (171) is positive definite, i.e., when 
|0) is a stable vacuum [l8|. H b can then be rewritten as 

H b = (H) + E " a h'lb' a + - A") , (9) 

a 

where A" stands for A* , u> a are the symplectic eigenvalues 
of T~i, i.e., the positive eigenvalues of the matrix 

^=( A A a+ -A?A + ).*<=(i-0 

(10) 

whose eigenvalues come in pairs of opposite sign (and 
which is diagonalizablc with real non-zero eigenvalues 
when T-L is positive definite), and b' a , are "collective" 
boson operators related to the local ones by a Bogoliubov 
transformation Z = WZ', i.e., 

with i^r) a and Qj) a the eigenvectors of M!H associated 
with the eigenvalues u a and — uj a respectively (such that 
W^MHW = MO,, with fl aa , = \uj a \5 aa ,). In order to 
preserve the boson commutation relations, which can be 
cast as ZZt - [(Zt) tr Z tr ] tr = M, W should satisfy 

WMW^ = M (12) 

which implies also W^MW = M and hence W^HW = 
fi. This entails WV-V tT U = 0, WU-V tT V = I, which 
are the natural orthogonality relations fulfilled by the 
eigenvectors of (|T0|) with normalization (^y a M.(^)ot = 1. 

The RPA matrix ((7J is of dimension 2n x 2n, with n 
the number of spins. RPA involves then an exponential 
reduction in the dimension (from (2s + 1)" to In for n 
identical spins). Moreover, in a translationally invariant 
system (sec. IIID|) . it can be further reduced to n 2 x 2 
matrices, becoming then fully analytic. 

B. The RPA ground state 



H b = (H)q + - V/ 'A - E A + b l b J + K A - 6 1 & ] + h - c -) The vacuum of the new bosons b' (b' a \0' b ) = 0) is 
= ( tf )o _i£ A *+i Z t HZ; (6) \0' b ) = C b exp[^Z ij blb]]\0 b ), Z = VU~\ (13) 

id 

Z=(t\, H= ( A A _A ^ \ , (7) where C b = (0 b \0' b ) = Dct^]" 1 / 2 is a normalization fac- 

V / V _ + / tor and Z a symmetric matrix. The associated RPA spin 

A% = ±y/s~s~ [J rxjx ± - i(J iyjx =p J ixjy )}, (8) state can then be defined as 

yij 

where Zt = (b\b) and A« = X i 6 ij . The choice of the |0 RPA ) = C s cxp[i ^ s i+Sj +]\0) . (14) 

mean field axes for the bosonization ((5]) ensures that no jjj 2^/sTsT 
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The expectation values generated by (fT4]) will be close to 
those obtained with the mapping ([5|), coinciding exactly 
up to second order in V (Appendix In contrast with 
|0), the state is entangled (unless V ^ 0). 

Let us note that for the quadratic Hamiltonian ([1]): 
i) |0rpa) = |0) if and only if |0) is an exact eigenstate of 
H, since H b contains the exact matrix elements connect- 
ing |0) with the rest of the Hilbert space: 

H\0) = {H) \0)-tJ2A i l\l i l j ), (15) 

id 

where |lilj) = ^±==|0) and we have used the mean 

field condition (li|if|0) = (h\hi\O t ) = (Eqs. ©-©)• 
Thus, if |0 RPA ) = |0), Z = and hence V = in W, 
implying A_ = 0. |0) is then an exact eigenstate by 
Eq. (fT5"| . Conversely, if |0) is an exact eigenstate, it is a 
solution of the mean field equations leading to A_ = 0, 
implying |0rpa) = |0) (although A + may be non-zero 
and / A"). In particular, when H has an exactly 
separab le g round state |0) (i.e., at the factorizing field 
BHHl), JOrpa) = |0). 

ii) |0rpa) is always exact for sufficiently strong fields 
(\B\ ^> J). In this limit |0) is the state with all spins 
Si fully aligned along — B l plus small corrections (A 1 w 
B l + sJ ■ B l /\B\ l ). Up to first order in A±/A, Eqs. 

(fj"0" ]) -([Tg |) lead then to Z ij w V l3 « j^r, entailing 

\0 RPA )^\0) + J2j^y\Ulj), (16) 

i<j 1 3 

which, by Eq. (|15[) . is just the first order expansion (in 
A_/A) of the exact ground state. 

In the case of a symmetry-breaking mean field, the 
RPA spin state allows to implement the necessary ro- 
tations for symmetry restoration: The exact ground sate 
will actually be close to the superposition with the correct 
symmetry of the degenerate RPA ground states (rather 
than to a particular RPA state), as will be discussed in 
Sec. Ill El in the context of parity-breaking. This restora- 
tion enlarges considerably the capabilities of the RPA. 

C. Bosonic evaluation of subsystem entropy and 
negativity 

The direct evaluation of many-body correlations and 
entanglement measures from the RPA spin state (fT4|) is 
in general difficult. However, the values of these quan- 
tities in the associated bosonic vacuum (|13p . which will 
be close to those obtained from (fT4|) . can be straight- 
forwardly evaluated using the general gaussian state for- 
malism [25I [26l | . The reduced density matrix of any sub- 
system is just a gaussian state, i.e., a canonical ther- 
mal state of an effective quadratic bosonic Hamiltonian, 
since Wick's theorem holds for the evaluation of the mean 
value of any observable, and in particular those of the 



subsystem. We may then evaluate its entropy and other 
invariants through standard expressions for independent 
boson systems. 

Let us formalize the previous scheme. We will use a 
generalized contraction matrix formalism, equivalent to 
that based on covariance matrices [25l |26| , which is more 
natural for the present RPA formulation. In the new 
vacuum |0' 5 ), (&„&„')<)' = ( b ' a b ' a <)o> = 0, implying 

Fij = (b]bi} , = {VV%, (17a) 

Gij = (&AV = (V^h ■ (17b) 

Eqs. ([5|)- (fT7|) determine the basic RPA spin averages and 
correlations, i.e., (s ifi } > = S^Fa - s») and, for i ^ j, 

(si+8j-)o> = 2y/siSj Fji , (si_s 3 -_)o' = 2y/SiSj G jt , 

(18) 

with (si±Sj Z )o> = 0, which coincide exactly with the av- 
erages derived from (|14|) up to second order in V, i.e., 
first order in the average occupation VV^ (normally very 
small outside critical regions). Through the use of Wick's 
theorem, we also obtain (si Z Sj Z )o' = (siz)o'( s jz)o' + 
\F l0 \ 2 + \G l] \ 2 iori^j. 
We may now define the generalized contraction matrix 

V={ZZ*)u-M= , (19) 

which exhibits the correct transformation rule under Bo- 
goliubov transformations: If Z = WZ' , then 

V = WD'W^ . (20) 

with V = (Z'Z'^o' -M. Eq. dTTJ) can in fact be written 
in the form (|20|) if W is the diagonalizing Bogoliubov ma- 
trix {IT]) and V the vacuum density (F' = G' = 0). We 
may then also obtain W and V through the symplcctic 
diagonalization of T>, i.e., through the diagonalization of 

VM = ( % _f_ p ) (21) 

such that W^VMW = VM, with V diagonal. 

Let us consider now a subsystem A of m < n sites. It 
will be characterized by a truncated contraction matrix 

V A = (Z A Z\) , ~M A =(l A Al ° + a Pa ) (22) 

where Z A contains just the bosons of sites in A. A sym- 
plectic diagonalization of V A will lead to 

V A = W A V' A W\, -D' A =( f Q 1 l u ) , (23) 

where f%*' = f^S aa ' with f% = VlX*)* ^ 
iV A M. A has eigenvalues f A and —1 — f and 
W A M A W\ = M A , with Z A = W A Z' A . The entangle- 
ment between A and its complement A is then given by 
the associated bosonic entropy, 

S{p b A ) = -T?p h A \og 2 p h A (24) 
= -Y,fl l0 S2 ft - (1 + fl) log 2 (l + f%) (25) 
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Here p\ = Tr^lO^O^ is the bosonic reduced density of 
subsystem A, which can be explicitly written as 

P \ = Cex v [-\Z\U A Z A ] = CcM-J2 LU A b 'lM^) 

ot 

where C = Yl a (^ + fa) anc ^ Ha, T^a are related by 

V A M A = [cMM A n A ) - 1]' 1 ■ (27) 

Here H. A represents an effective "Hamiltonian" matrix for 
subsystem A with symplcctic eigenvalues u} A such that 
f° = (e w ° - l)- 1 (and hence -1 - f% = (e~ u ^ _ 
Eq. (f2"6")) leads then to the contraction matrix (|22[) . and 
hence to the same expectation values as the full vacuum 
f° r an y operator of subsystem A. 

Eq. |25|) provides a tractable RPA estimation of the 
entanglement entropy of any subsystem. It is shown in 
the Appendix [X] that a direct spin evaluation of the sub- 
system entropy based on the RPA state (Tl4|) coincides 
with (|2"5|) up to second order in V. 

On the other hand, the internal entanglement of sub- 
system A with respect to a partition (B, C) of A (where 
the complement A plays the role of an environment) can 
be measured through the corresponding negativity [l3{, 
defined as minus the sum of the negative eigenvalues of 
the partial transpose p l A of p A : 

N BC = %(Tr\p t ?\-l). (28) 

Expectation values with respect to (p b A ) tc of an observ- 
able O h A correspond to those of the partial transpose 
(0 A ) tc with respect to p\. This implies the replace- 
ments Fij <-> dj, Fj'j <-> Fjj', Gj'j <-> Gj'j, in the con- 
traction matrix for j,j' G C, i G B, leading to a matrix 
T) A with symplcctic eigenvalues f A . The latter can now 
be negative. We may then still write {p h A ) tc as in Eq. 
(|26p in terms of an effective matrix r hi A with symplectic 
eigenvalues lo a such that f A = {e" A — 

Since the trace remains unchanged (Tr {p b A ) tc = 1), 
\e-**\ < 1, implying f% > -1/2. A negative f% > -1/2 
corresponds to e~^ A < and hence to a non-positive 
(p h A ) tc , indicating an entangled p b A with respect to this 
bipartition. We then obtain, noting that (1 + e~" A ) _1 = 
(1 + f%)/{l + 2f%), the final result [il,^,^ 

n<o 1 + 2fA 

which allows the evaluation of the negativity (|28p . Nega- 
tivities obtained from the spin density matrices coincide 
with this result up to first order in V (Appendix [A"|l . 

In the case of a global bipartition (A, A), N A j b ecomes 
a function of the reduced density p A , namely jl7| 

N AA = |(Tr||0>(0|^| - 1) = ±[(Tr^) 2 - 1] . (30) 

In a boson system, this implies that N AA , a limit case of 
Eqs. (|28p- (|2TJ)) . can be also expressed just in terms of the 



symplectic eigenvalues f A of the contraction matrix D A : 

N AA = k[\{(Vn + ^+n) 2 -^- (3i) 



D. Translationally invariant systems 

The only quantities required in the bosonic RPA 
scheme are, therefore, the basic contractions (fTTjl . Their 
evaluation becomes remarkably simple in translationally 
invariant systems, either in one or d dimensions, i.e., 
systems with a common spin Sj = s in a uniform field 
B l = B with couplings dependent just on separation: 

J^ ju = J** v {i-j) (32) 

where J» V {1) = J vtl (-l), and J^(-l) = J^(n - I) in 
a finite cyclic chain or system (in d dimensions, i,j,l,n 
stand for d-dimensional vectors). We will also assume a 
uniform mean field A 1 = A, which should then satisfy 

A" = B" J o V Mo > 4" = E J/ "« ' ( 33 ) 

V I 

with (s)o = — sA/A (Eq. ([4])). The uniform mean field is 
thus determined just by the total strengths Jq V . 

Choosing again the z axis in the direction of A, such 
that (s i(ti ) = -sS^ and B^ + sj£ z = W, with A > 0, 
the bosonized Hamiltonian will have the form ^ with 
couplings A !j? = A± (i—j)- By means of a discrete Fourier 
transform of the boson operators, we can rewrite it as 

H b = (H) + £(A - A k + )b\b k ±(A k _b\bl k + h.m 
k 

n-1 

A k ± = ^V 27rfcZ/n A ± (0, (35) 

1=0 

where k = 0, . . . , n - 1 and b k = fa £" =1 ^ ki/n bj are 
boson operators in momentum space, with b-k = b n - k . 
Diagonalization of (|34|) is straightforward and leads to 

H b = (H) + E co k b'lb' k + \{u k A + A* ) , (36) 
fe 

where w fe = Co k - \{A\ - A^ fc ), b'\ = u k b\ + v k b^ k and 

u k = ^(X-A k ir-\A k ^ , (37) 

/ A - A*. + u k A k _ l\-A k -Co k , 

Uk = V 2^ ' Vk = |#| V 2^ (38) 

with A^ = i(A^ + A; fe ), u\ - \v k \ 2 = 1, and u k = u_ k , 
v k = v^ k . All uj k should be real and positive for a stable 
mean field, implying the stability conditions 

< |A*| < A- A*, fc = 0,...,n-l. (39) 
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We can now obtain the basic contractions explicitly, 

t A - 
(&fc&fc')o' = $kk'\v k \ , (bkb-k')o' = $kk'U k v k = ^fc( 4 °) 

which lead finally to (Eq. (fT7]l) 

Fa = F(i -3) = lE e- a ^ )/n \v\\ , (41a) 

k 

Gij = G(i -j) = -V e-' ;27rfc ( l -^/" Ufc ^ . (41b) 
fe 

For strong fields |£?| such that A > |A±|, Wfc^fc sa ^A^/A 
and |v|| w i|A*L| 2 /A 2 . The RPA vacuum ([13]) becomes 

|0' b ) = C h exp[± £ Z(i - j)blb]]\O b ) , (42) 

i,j 

where C 6 = IL u~ 1/2 and Z(i) = £ £ fe e^ 2 ^™^. 

Thus, these systems allow an analytic evaluation of the 
contractions (fl"7|) . Both the mean field equations (j33|) 
and the RPA Hamiltonian (|34|) become independent of 
the common spin s after a rescaling J^ u (l) — > J^(l)/s, 
which we will adopt in what follows and which indicates 
that RPA is describing the large spin limit of the system, 
as is apparent from Eq. (|5|). 



E. XYZ systems 

Let us now examine in more detail the previous formal- 
ism in a translationally invariant spin s array with XYZ 
couplings of arbitrary range in a uniform transverse field: 

H = B s lz - ^ J ^ i ~ -?') s v s j> ■ ( 43 ) 

i i^j fi=x,y,z 

Eq. (|43[) commutes with the S z spin parity, 
[H,P z )=0, P, = exp[*7r 

i 

for any value of its parameters, such that the exact 
ground state in a finite array will always have a defi- 
nite parity outside degeneracy points. We will focus here 
on the ferromagnetic type case where J x (l) > V I with 

\Jy(l)\ < Ml) , (44) 

which exhibits a normal and parity breaking phase at the 
mean field level. 



1. RPA around the normal state 

For the Hamiltonian (|4"3")l . the state |0) with all spins 
fully aligned along the — z axis is always a solution of the 



mean field equation (|33|) . being the lowest solution for a 
sufficiently strong field B. It leads to A 1 = \5^ z , with 

A=|S| + J z °>0, J° z =J2Ul)- (45) 

i 

All previous equations can then be directly applied. Now 
A± ( Z ) = J*d)±Jy(i) = A±(-Z), implying A k ± = A ± k and 

co k = y^A - J k ){\ — Jy) , (46) 

where J k = e a * kl l n J M (/) {A k ± = ^^-). This solu- 
tion is therefore stable provided J k < A Vfc and [i = x,y, 
i.e. for \B\ above a certain critical field B c . In the case 
(|4"4f . the strongest condition is obtained for k = 0, i.e., 

\B\ > B c = J° x - J° . (47) 



2. RPA around the parity breaking state 

For \B\ < B c , the normal state becomes unstable: the 
lowest normal RPA frequency w° vanishes for \B\ — > B c 
and becomes imaginary for \B\ < B c . The lowest mean 
field for \B\ < B c corresponds instead to a parity- 
breaking state with all spins aligned along an axis in the 
xz plane forming an angle 9 with the z axis: 

|0) -> |9) = |0i . . . 6 n ) , \6j) = expH^HO,-) . (48) 

This leads to (sj)o = — s(sin 9, 0, cos 9) = — sA/A, with 

A = J°, cos9 = B/B c , (49) 

as determined by (|33|) . We should now express the origi- 
nal spin operators in terms of the rotated operators, i.e., 

Six = s ix f cos 9 + s lz > sin 9, s iz = s iz , cos 9 - s ix , sin 9 

(50) 

with Si y = Siyi. The RPA around this state amounts 
therefore to the replacements 

A J°, J k -> J' k x = J k cos 2 9 + J k sin 2 9 , (51) 

in Eq. with J k unchanged and A^ = \{J' k x ± J k ). 

Correlations (s^'Sj^rpa of rotated spin operators 
have the same previous expressions (|17j) , whereas those of 
the original operators must be obtained using Eqs. (|50[) . 
It should be remarked, however, that in a finite system, 
the associated RPA spin state will no longer be a good 
approximation to the actual ground state due to par- 
ity breaking. Parity restoration, at least approximately, 
must be implemented before obtaining final results. We 
will not discuss here the case of a continuous broken sym- 
metry (arising for instance in the XX Z case), which can 
bc treated through the RPA formalism of ref. @. 
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3. Definite Parity RPA ground states 

Since [if, P z ] = 0, the parity breaking mean field state 
|0) is degenerate: Both |0) and | — 0) = P z |0) are mean 
field ground states. In order to describe the definite par- 
ity ground states, the correct RPA ground state should 
be taken as the definite parity combinations 
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10 



RPA, 



±1-0 



RPA/ 



RPA/ 



y/2(l ± (-0rpa|0rpa)) ' 



(52) 



where | ± 0rpa) are the RPA states around each mean 
field. The overlap (-0rpa|©rpa) = (0rpa|-P z |0rpa) is 
proportional to the overlap between the two mean fields, 



0|0) = cos 2ns 9 = (B/B c ) 



2ns 



(53) 



which is small except for B — > B c or small ns. 

Neglecting the previous overlap, Eq. (|52|) will lead to 
reduced densities 



pi 



\[ PA {e) + PA {-6)\ 



(54) 



~©rpaI©rpa) 



oc 



provided the complementary overlap 

(lf~) 2 ' n_ "' 4 ^ s can also b° neglected. Here ^(±0) are 
the reduced spin densities determined by each RPA state, 
given up to 0(V 2 ) by the expressions of Appendix IaI 

The restoration (|54|) is essential to achieve a good 
description of the actual subsystem entropy, although 
its main effect for a not too small subsystem A is ac- 
tually quite simple: If the product P a(@)pa(—®) oc 
(B/B c ) 2nAS can be neglected, Eq. (|54|) can be consid- 
ered as the sum of two densities with orthogonal support 
and identical distributions, leading to 



S(P A ) 



S( PA (6)) + 1, 



(55) 



where S(p A (Q)) can be evaluated through the boson ap- 
proximation (I25p . Under the same assumptions, the ef- 
fect on the global negativity ([30)) is just 



N A a(Pa) « 2N AA ( PA (0)) + 



(56) 



as Try p A « v^Tr^/ p A {9), while the subsystem nega- 
tivity Nbc °f a bipartition (B, C) of A remains approx- 
imately unchanged: N BC (p A ) « N B c(pa(9))- 

When the product P a(Q)pa(~ 0) cannot be neglected 
(as in a subsystem of two spins), we should in principle 
construct the spin density (|54[) . This can be done by 
rotating P a {9) (Eq. (|A2p in the mean field frame) to the 
original z axis and removing all parity breaking elements 
(which is the final effect of Eq. For instance, the 

reduced two-spin density for s = 1/2 has the blocked 
form (|A2[) in the standard basis of Si Z Sj Z eigenstates in 
the normal phase as well as in the parity breaking phase 
after parity restoration [l2j]. The final effect on S(p A ) is 
the replacement of the term +1 in (f55"j) by the entropy of 
the reduced mean field mixture — Ylv=± Qv 1°S2 If , with 
q± = |(1 ± (B/B c ) 2sa ), plus small RPA corrections. 



While p^ are both identical in the approximation (|54[) , 
the actual p^ in a small system will depend on parity. 
The correct parity in such a case should be chosen as that 
leading to the lowest energy E^ pA = (0rpaI-^1©rpa)- 



4- Factorizing Field 



The explicit value of the basic RPA couplings A± in 
the parity breaking phase are, using Eqs. (f5"Tj) - (|49)) . 

A^ 



|[(J^-J*)(B/S C ) 2 + J*±J*)] (57) 



In the case of a common anisotropy, such that the ratio 

, _ Jy(l) - Ml) 



Mi) - Mi) 



(58) 



is independent of the separation I, we have Jy — J. 
X(J* - J k z ) and hence A^ = - J k z )[{B/B c ) 2 - 

It is then seen that if x € [0, 1], A* = V k when 



\B\ = B s = B c ^ 



(59) 



with all A* changing sign at \B\ = B s . Here B s is the 
factorizing field [2j, 1 12t |22h24| : At B = B s the parity 
breaking mean field state becomes an exact ground state, 
since the RPA corrections vanish (sec. Ill Bp . This effect 
is independent of the number of spins n (as long as \ is 
constant) and spin s (with the present scaling) . Nonethe- 
less, the actual side limits at B = B s will be given by 
the definite parity states (|52|) . which are still entangled. 
As a consequence, the subsystem entropy S(p A ) and the 
negativity N AA will actually approach a finite value for 
B — > B s (1 and 1/2 respectively in the approximation 
(|55|) - (|56j) ). while the entanglement between two spins 
will reach there infinite range fiol - fl^j ]. Note finally that 
at B = B s , A^ = Jy and hence, 

u k = J° - J* . (60) 



III. APPLICATION 



A. Spin s pair 



As a first example, let us consider a system of two 
spins s coupled through the Hamiltonian (|4"3"| . We can 
obviously always set here J x > \J y \ (Eq. (|44|) ), since the 
sign of J x can be changed by a 7r-rotation around the z 
axis of one of the spins (and we can always set \ J X \ > \J y \ 
by a proper choice of axes). The Fourier transform of 
J/j,(l) = <5/iJ M reduces here to = ( — l) fc J M , k = 0,1, 
leading to an attractive and a repulsive normal mode: 



LUQ = J (A - J X )(X - Jy) , UJl = J{X + J X )(X + Jy) 



The contractions (|4T|) become Fu = 



2<5 ij) 2^y'' 



4w 



A-A + A+A + 
4u) 



pa 



— (1 — 25ij), where A± 
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FIG. 1. Entanglement between two spins s as a function of the 
transverse field B for an XY coupling with J y /J x = 0.5. The 
exact entanglement entropy Se = S(pi) (top) and negativity 
(bottom) for different values of the spin s, and the bosonic 
RPA results, Eqs. (|61[) . ()63|) are depicted. The exact results 
approach those of RPA as s increases, differences for not too 
small s arising just for B close to B c = J x . At the factorizing 
field B s w B c /V2, Ss = 1 while A 1 ^ = 1/2. 
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FIG. 2. Top: Left: The entanglement entropy obtained from 
the definite parity RPA spin state (|52|) (dashed-dotted line), 
compared with the bosonic RPA result (|61[) and the exact 
value, for s = 10 at the same parameters of fig. [1] The result 
from the RPA spin state improves the bosonic RPA for B just 
below B c - Right: The average local boson occupation /, Eq. 
(162[) , which is small away from B c , and the negative eigenvalue 
/ of the partial transpose of the contraction matrix (/ ps \/J 
for small /). Bottom: Left: RPA energies tJo, ui, together 
with the mean field energy A and the mean RPA energy ZJ 
appearing in (|62[) . Right: The quantities Zk = Vk/ut for 
k = 0, 1, which determine the RPA state (|42[) and vanish at 
the factorizing field B s . 



\{J X ± J y ) and replacements ((51")) are to be applied for 
\B\ < B c . The ensuing entanglement entropy of the pair 
in the bosonic approximation (|25| is just 

S( Pl ) = -/log a / + (1 + /) log 2 (l + /) + 6, (61) 

/=2 V 1 + 1 )> U = o 62 

where / = J (Fu + i) 2 — (Gn) 2 — \ is the positive sym- 

plectic eigenvalue of the 2x2 contraction matrix for one 
spin and S = (1) for \B\ > B c (< B c ) in the approx- 
imation (|55| . valid for (B/B c ) 2s <C 1). For small /, we 
may just use S(pi) w /(log 2 e - log 2 /), with / w Fn, in 
agreement with the results of Appendix A. 

Thus, at the RPA level entanglement is determined by 
the average local occupation / and driven by the ratio 
- — — , which is small away from B c and vanishes at B = 
B s (where ZJ = A = J° by Eq. (J60J), and hence / = 0). 
For \B\ > B c , f {^j^-f, while in the vicinity of 

B„ f cx (fl - B s ) 2 . For H- B c , / « |^/^ « 

iS-Sel- 1 /^ W ith 5(pi)pslog 2 /e. 
The bosonic RPA negativity p8 |) -([29 |) becomes 



N12 = =f+ Vfif + 1) (63) 

where / = / — ^/ /(/ + 1) is the negative symplectic 
eigenvalue of the 4x4 contraction matrix. Correction 



1(55)1 (A 2 i 2A 2 i + |) should be applied for |S| < B c . 
For small /, we have simply A/i 2 ps — / ps ^/J. This will 
lead to a slope discontinuity of Ai 2 at the factorizing 
field B s (see Fig. QJ, as / vanishes there quadratically 
(Ni 2 - | oc \B - B s \ for B ps B g ). On the other hand, 
for / -> oo -> S c ), /-> -|, with / ps -I + ^ and 
Ai 2 ps 2/. Both S(pi) and iVi 2 are concave increasing 
functions of / and measure the entanglement of the pair. 

Comparison with exact numerical results, obtained 
through the diagonalization of H (a (2s + l) 2 x (2s + l) 2 
matrix), are shown in Fig. [1] for the XY case (J z = 0) 
with anisotropy \ = Jy/Jx = 0.5. Exact results are seen 
to rapidly approach the RPA values (f6"l"]) -(f63 j) as the spin 
s increases, the discrepancy for finite s arising just in the 
vicinity of B c or for very small s, i.e., where tunneling 
effects arising from the non-zero overlap ((53)) between the 
degenerate parity breaking states become appreciable. 

Nonetheless, this overlap can be taken into account 
using the full definite parity RPA spin state ((52)) with 
lowest energy, which for finite s improves results for B 
close to B c (but otherwise yields results almost coincident 
with those of the corrected bosonic RPA), as seen in fig. (2) 
Eq. (l52l) also yields the exact side limits at the factorizing 
field [l2|] for any s, although for x = 0.5 these limits 
rapidly approach the high spin values S(pi) = 1 and 
N12 = \ predicted by the approximations ((55)) - ((56)) . 

Fig. 2 also depicts the behavior of the average occu- 
pations / and /. The former is seen to bc quite small 
(/ ^ 0.05) except in the vicinity of B Cl implying that 



away from B c , all bosonic RPA results can be reproduced 
by the spin densities of Appendix A, with / w y/J. In 
the bottom panels we depict the RPA energies ljq, uj\ and 
the RPA state coefficients Zk = Ufc/ufe used in Eq. (|4"2"j) . 
Although luq vanishes at B c , the difference A — uJ, respon- 
sible for entanglement, remains everywhere quite small. 
Both Zk vanish and change sign at the factorizing field 
B s , indicating a qualitative change in the type of correla- 
tions at this point: Entanglement between two spins 1/2 
is well known to change from antiparallel to parallel (in 
the original frame) at B s [l2j , an effect arising within the 
RPA from this sign change. 



B. Fully connected spin system 

Let us now consider a fully and uniformly connected 
XYZ array of n spins, where 



Ml) = (1 - *io)V(n - 1) ; 



(64) 



in (|43j) . This scaling ensures a finite intensive energy 
(H) jn for large n and finite J M . Entanglement prop- 
erties of this well-known model [HI, HtJ for s = 1/2 in 
the large n limit have been previously analyzed [281 ] , in- 
cluding recently Holstein-Primakoff based bosonization 
HE El E3]. Direct application of the present RPA 
formalism will be here shown to yield full analytic ex- 
pressions for any size n and spin s. The present treat- 
ment docs not exactly coincide with that of refs. fl6L f2~T| . 
since the absence of self-interacting terms oc s^Sj^ (non- 
trivial for s > 1/2) is here exactly taken into account and 
leads to repulsive RPA corrections (cJi), non-zero for fi- 



nite n. The Fourier transform of 



J, 



-J M /(n — 1) for k = 1, 



I64J) 

- 1. 



is J° = 



J,, and 



leading again to 



two distinct RPA energies: One associated with a funda- 
mental attractive mode (loq) and n — 1 degenerate weak 
repulsive modes Uk = u>i, k ^ 0, which just add a small 
repulsive correction accounting for the absence of self- 
energy terms: 



w 



y/(X-J x )(X-J y ), Wi = 

where the replacements (foTj) are to be used for B < B c . 
The ensuing contractions (|41|) become here obviously in- 
dependent of separation for i ^ j: 



r% 3 ~ 2n I O) 



A — A; 



(1 - nSij)] 



±<5- 

2 U »3 ' 



Gn = -±-[—-—(l-n6ij)], 



(65a) 
(65b) 



and imply that for any bipartition (L, n — L), the entan- 
glement entropy S(pl) will depend just on L. Moreover, 
there is again a single non-zero eigenvalue Jl of the re- 
duced matrix T>l of L spins for any L (see Appendix B), 
such that in the bosonic approximation (|25|)- (|55l). 



S(pl) = -h log 2 Sl + (1 + h) log 2 (l + h) + S , (66) 
f L = i[ v / l + 2a L A - 1] , a L = L{n - L)/n 2 (,67) 
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FIG. 3. Results for the fully connected spin 1/2 array of 
n = 100 spins. Top: Left: Exact entanglement entropies 
Se(L) = S(pl) of subsystems with L < n/2 spins as a func- 
tion of the magnetic field. Right: Comparison between exact 
and RPA results for S(pl)- Bottom: Left: Exact and RPA 
results for S(pl) as a function of the subsystem size L at four 
different field ratios B/B c . Right: Magnetic behavior of the 
average boson occupation number Q67p f° r L — 25 and the 
negative symplectic eigenvalue (I70[) of the partial transposed 
contraction matrix for different L, m. 



where 5 = 0(1) for \B\ < B c {{B/B c ) 2Ls < 1) and 
n 2 (A 2 — uJ 2 ) _ uj + (n — l)wi 



A 



2(n — l)woWi 



(68) 



For n = 2 we recover Eqs. (|61|) (|62l) . while for large n, 

A-A° 

A « — — 1. Entanglement is then driven again by the 
\ 2 — 2 

ratio ~" , which is small away from B r and vanishes at 

B s . For small A, / L a ±a L A, with A « f Qj^ ^ ) 2 

for |B| > B c and A oc (B-B s ) 2 in the vicinity of B s . For 
B -> B c , f L oc ^(B - Be)" 1 /* and S(p £ ) « log 2 / £ e. 

The bosonic negativity of a bipartition (771, L — m) of 
a subsystem of L < n spins can again be explicitly ob- 
tained, since there is also a single negative eigenvalue 
fhm of the partial transpose of the contraction matrix 
(see appendix |B|) : 



N„ 



,L-r 



h 



Lm 



1 + 2/. 



Lm 



1 + 7Lm A - y/8p Lm A + Y Lm A- 
f 4/3z, m , /3i m = m(L - m)/n 2 . 



7im = «i "I" *Pim 

For a global partition (L = n), a T 



(69) 

(fo) 

(71) 



while j3„ 



and /„£ = f L - y/f L {f L + 1), with = / £ 



\J + 1), as in Eq. (f6"3) . In general, for small A, 



such that for strong fields, /l^ 



(72) 



a/^t^t^, while 



l n-l 4B 

for B close to B;,, /i m oc \fPi~^\B — B s \. On the other 



9 






111)111 

Exact 

CPA 




i 


. . . 1 . . . 



ill; 


— I — I — I — ■ — i — i — i — | — i — i — i — 

\\\ Exact 




;/ \ RPA \ 

ij \ 




ij \ 

UA \m=10 










. . . 1 



0.5 1.0 1.5 B/J> 

FIG. 4. Top: Left: Exact global negativities N(L) = N L ,n-L 
between L and n — L spins in the fully connected array. Right: 
Comparison between exact and RPA results for N(L) for two 
values of L. Bottom: Left: Exact subsystem negativities 
Nm t L~m between m and L — m spins in a subsystem of L = 20 
spins. Right: Comparison between exact and RPA results for 

m,L — m- 



hand, for B — S- B c , f 



Lm 



0(\B 



a L +4,p L , 

Bel 1 / 2 ) if ax 7^ 0, implying that subsystem negativities 
N m .L-m with L < n remain finite at B c (in agreement 
with the results of [IH), as fi,m remains above — |. 

In the parity breaking phase, the replacement (|56p 
(N — > 2N + |) should be used for global negativities 
N„.L~n, whereas subsystem negativities N m ^- m remain 
unchanged after parity restoration if (B / ' B c ) 2s ^ n ~ L ^ and 
(B/B c ) 2sL can both be neglected. 

Eqs. (|66|) - ([70|) represent essentially the exact expres- 
sions for the subsystem entropy and negativity for large 
spin at finite n, as well as for large n at finite spin, as ver- 
ified by exact numerical calculations. For instance, exact 
(obtained through diagonalization of H) and RPA results 
for a spin 1/2 array of n = 100 spins are shown in figs. 
[3] |U RPA results for the entanglement entropy are quite 
accurate except in the vicinity of B c , differences decreas- 
ing as n or s increases. For large L they were obtained 
with the previous expression (|66|1 whereas for small L 
(like the L = 2 case) , we have used the proper spin state 
(|54|) . whose main effect is to take into account the correct 
overlap for B below but close to B c (roughly, S replaced 
by the entropy of the reduced mean field superposition). 

The variation of S(pl) with L at fixed field (bottom 
left panel in Fig. [3} is also correctly predicted, being 
quite accurate both in the normal and parity break- 
ing phase for fields not too close to B c . The bot- 
tom right panel shows that fi, remains small except for 
B around B Cl while jhm becomes also small as L de- 
creases, in full agreement with Eq. ([72jl . RPA results for 
global (N nt L- n ) and in particular subsystem negativities 
{N m ,L-m for L < n), which are much smaller and vanish 
at B s , are also very accurate, as seen in Fig. [5] Sub- 
system negativities were directly obtained with Eq. (|69p , 
whereas global negativities were corrected with Eq. (|56p 



for B < B c and large L and using ([54]) for L = 2. 



IV. DISCUSSION 

We have shown that the mean field plus RPA method 
is able to provide, through the bosonic representation, a 
general tractable method for estimating, in the ground 
state of general spin arrays, the entanglement entropy of 
any bipartition of the whole system as well as the negativ- 
ity associated with any bipartition of any subsystem. The 
approach becomes fully analytic in systems with transla- 
tional invariance, where no numerical diagonalization is 
required for obtaining the basic contraction matrices. 

The bosonic treatment provides essentially the exact 
behavior of the system in the large spin limit. Finite 
spin corrections can be taken into account through the 
corresponding RPA spin state, which allows in particu- 
lar to implement the non-negligible symmetry restoration 
effects in the case of the parity-breaking mean field, but 
which otherwise yields result which are in full agreement 
with the bosonic treatment at first order in the average 
local boson occupation. The latter is normally very low 
away from critical regions. 

Through direct application of the present method, sim- 
ple analytic expressions for the entanglement entropy and 
negativities for a spin s pair and for a fully connected 
array of n spins s in a uniform field, have been straight- 
forwardly obtained, which depend explicitly on the RPA 
energies. The agreement with exact numerical results 
is confirmed to improve as the spin s increases at fixed 
size, and in the fully connected case also as n increases 
at fixed s, differences being in fact negligible away from 
the critical region for not too small s or size. 

An important general prediction that arises from the 
present treatment is that entanglement from elementary 
excitations approaches a non- vanishing spin independent 
limit as the spin increases. An RPA quantum regime, 
characterized by weak entanglement, emerges then be- 
tween strictly classical and strongly quantum regimes. 

The authors acknowledge support from CIC (RR) and 
CONICET (JMM,NC) of Argentina. 



Appendix A: RPA spin densities 

We will here construct the spin density matrices com- 
patible with the RPA spin state JT3J and the contractions 
(|17p up to second order in V, i.e., first order in the av- 
erage occupation VV^ (implying zero or one boson per 
site). At this order, F GG t (Eqs. ([Tfjl l and the sup- 
port of p = |0rpa)(0rpa| is just the subspace spanned 
by the mean field state |0) plus the two site excitations 
|l<lj) (Eq- OS)), leading to 



Gt 1 



G 

G^G 



(Al) 
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where G denotes a column matrix of elements Gy , i < 
j. At this order, p 2 = p. The ensuing reduced density 
matrix pa = Tr^ p of a subsystem A of L spins becomes 

( G A G\ G A \ 

PA~\ F A - G A G\ (|2) 

\ G\ 1-trFA + G^GU/ 

where Fa, G A are the L x L contraction matrices of 
subsystem A and G A the concomitant column vector 
(of length L(L — l)/2). The central block contains the 
one-site elements |li)(lj| arising from the partial trace 
of GGL Here we have used the approximate identity 

Efc G A GikG{j « Fij-T,keA G lk G jk for », J e A (and ne- 
glected diagonal elements Gu , of higher order due to the 
absence of self-energy terms), which allows to write p A 
entirely in terms of local contractions. Eq. (|A2|) is then in 
agreement with direct state tomography at this order (for 
i,j,k,l € A, - blb k )) , « (F A - G A G A ) 13 , 

(6j6]6 ft 6,) , w GkiGij). Up to 0(U 2 ), pa is a positive 
matrix with Trp A = 1, but is no longer pure. 

Its entropy S(p A ) = — Tip A log 2 p A is determined, at 
this order, by the central block p\ = F A — G A G A , 

S{p A ) « tr p\ (log 2 e - log 2 p\) , (A3) 

which coincides with Eq. (|25[l up to second order in V 
(at this order f A coincides with the eigenvalues of p\ and 
Eq. m becomes « £ Q /^(log 2 e - log 2 /£)). 

On the other hand, the leading term in the negativ- 
ity arising from a bipartition (B, C) of A is of first order 
in V and is just the sum of the singular values of the 
submatrix Gbc (of elements Gy, i £ B, j G G), whence 
-ZVbc ~ tr [GscGjgp] 1 / 2 . At this order, the negative sym- 
plectic eigenvalues f A in (|29j) arc again minus the singular 
values of G_bc, while Eq. (|28|) becomes Asc ~ — S a /ji, 
leading again to the previous result. 

Let us finally notice that Eq. (|A2j) always commutes 
with the S z parity (along the mean field axis) of subsys- 
tem A, i.e., [p A ,P zA ] = 0, P zA = exp[iTTj2 ieA (siz - s»)]. 
In the case of two spins G A has length 1 and Eq. 
(|A2j) is just a 4 x 4 blocked matrix, while in the case of 
a single spin i, G A has length and Eq. (|A2[) becomes 
just Pl a F«|li>(li| + (1 - JiOlOiXOi). 



Appendix B: Fully connected system 

In the fully connected XYZ spin system, the con- 
tractions (pj|) are of the form = F dy- + F 1; Gjj = 
G <5 i:( - + Gi, with Fo,Fx, G ,Gi real. The ensuing con- 
traction matrix T>l for a subsystem of L spins will then 
have symplectic eigenvalues (see also [26| ) 

/l = ^(Fo + LF + |) 2 - (Go + LGi) 2 - i (Bl) 

/o = v/(F + i) 2 - G 2 - | (B2) 
plus their partners 1 + /l, 1 + /o, where Jl is non- 
degenerate while fo has L— 1 degeneracy. Eqs. (|B1|) (|B2|) 
can be obtained either by a Fourier transform of the lo- 
cal operators or by noticing that the L x L contraction 
matrix Fl can be written as Fl = Fo/l + FlI^I^ (and 
similarly for Gl), with Tj, the L x L identity and 1l a 
column L x 1 vector with unit elements. Fl and Gl will 
then be diagonal in the same local basis with eigenval- 
ues Fo + LF\ and Fo (L — 1 degenerate) and similarly 
for Gl, which leads to Eqs. (|B1|) (|B2[) . In the case of a 
global vacuum, fo — (since for L = n, we should have 
fh=n = fo = 0), implying a single positive eigenvalue /l 
for any L < n. Eq. (|BT|) leads then to Eq. (|67|). 

For evaluating the negativity N mp of a bipartition 
(m,p) of a subsystem of L spins (m + p = L), we 
may first note that Fl will be composed of blocks 
F — FnT 4- F 1 1* F — Fl 1* = F t and 

^ mm — 1 [) J m T i J- m , mp — ± 1 x m x p — x pm <aA1,a 

F pp = FqI p + Filplp, and similarly for Gl- A local 
transformation allows then to write Fl as a direct sum 
of a (L — 2) x [L — 2) diagonal block Fo/l-2 plus the 
block F / 2 + Fi(^j^), and similarly for G L . The en- 
suing partially transposed contraction matrix will then 
have symplectic eigenvalues fo = fo (Eq. (|B2[0 . L — 2 
degenerate (with f = for a global vacuum) and 

f~L = l\/ T ^ 2 ± V(Tr^l 2 ) 2 - 16dcl4 - \ (B3) 

together with their partners 1 + /o, 1 + /i m , where 
^ = <i FG ~j GF ) is a 4 x 4 matrix with blocks A FG = 

(| + F )/ 2 + QtkpgTpF! 1 ) and similarly for ^4gf- Here 
/lm > but f£ M < 0. The latter is the single negative 
symplectic eigenvalue given in Eq. (|70| . 
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